{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First Order Local Uncertainty Analysis for Chemical Reaction Systems\n",
    "\n",
    "This ipython notebook performs first order local uncertainty analysis for a chemical reaction system\n",
    "using a RMG-generated model.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from rmgpy.tools.uncertainty import Uncertainty\n",
    "from rmgpy.tools.canteraModel import getRMGSpeciesFromUserSpecies\n",
    "from rmgpy.species import Species\n",
    "from IPython.display import display, Image\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define the CHEMKIN and Dictionary file paths.  This is a reduced phenyldodecane (PDD) model.\n",
    "\n",
    "# Must use annotated chemkin file\n",
    "chemkinFile = 'uncertainty/chem_annotated.inp'\n",
    "dictFile = 'uncertainty/species_dictionary.txt'\n",
    "\n",
    "# Alternatively, unhighlight the following lines and comment out the lines above to use the minimal model,\n",
    "# which will not take as long to process\n",
    "# Make sure to also uncomment the specified lines two code blocks down which are related\n",
    "\n",
    "# chemkinFile = 'data/minimal_model/chem_annotated.inp'\n",
    "# dictFile = 'data/minimal_model/species_dictionary.txt'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the `Uncertainty` class object with the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uncertainty = Uncertainty(outputDirectory='uncertainty')\n",
    "uncertainty.loadModel(chemkinFile, dictFile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We can now perform stand-alone sensitivity analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Map the species to the objects within the Uncertainty class\n",
    "\n",
    "PDD = Species().fromSMILES(\"CCCCCCCCCCCCc1ccccc1\")\n",
    "C11ene=Species().fromSMILES(\"CCCCCCCCCC=C\")\n",
    "ETHBENZ=Species().fromSMILES(\"CCc1ccccc1\")\n",
    "mapping = getRMGSpeciesFromUserSpecies([PDD,C11ene,ETHBENZ], uncertainty.speciesList)\n",
    "\n",
    "initialMoleFractions = {mapping[PDD]: 1.0}\n",
    "T = (623,'K')\n",
    "P = (350,'bar')\n",
    "terminationTime = (72, 'h')\n",
    "sensitiveSpecies=[mapping[PDD], mapping[C11ene]]\n",
    "\n",
    "\n",
    "# If you used the minimal model, uncomment the following lines and comment out the lines above \n",
    "\n",
    "# ethane = Species().fromSMILES('CC')\n",
    "# C2H4 = Species().fromSMILES('C=C')\n",
    "# Ar = Species().fromSMILES('[Ar]')\n",
    "# mapping = getRMGSpeciesFromUserSpecies([ethane, C2H4, Ar], uncertainty.speciesList)\n",
    "\n",
    "# # Define the reaction conditions\n",
    "# initialMoleFractions = {mapping[ethane]: 1.0, mapping[Ar]:50.0}\n",
    "# T = (1300,'K')\n",
    "# P = (1,'atm')\n",
    "# terminationTime = (5e-4, 's')\n",
    "# sensitiveSpecies=[mapping[ethane], mapping[C2H4]]\n",
    "\n",
    "# Perform the sensitivity analysis\n",
    "uncertainty.sensitivityAnalysis(initialMoleFractions, sensitiveSpecies, T, P, terminationTime, number=5, fileformat='.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Show the sensitivity plots\n",
    "for species in sensitiveSpecies:\n",
    "    print '{}: Reaction Sensitivities'.format(species)\n",
    "    index = species.index\n",
    "    display(Image(filename=os.path.join(uncertainty.outputDirectory,'solver','sensitivity_1_SPC_{}_reactions.png'.format(index))))\n",
    "    \n",
    "    print '{}: Thermo Sensitivities'.format(species)\n",
    "    display(Image(filename=os.path.join(uncertainty.outputDirectory,'solver','sensitivity_1_SPC_{}_thermo.png'.format(index))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to run local uncertainty analysis, we must assign all the uncertainties using the `Uncertainty` class' `assignParameterUncertainties` function. `ThermoParameterUncertainty` and `KineticParameterUncertainty` classes may be customized and passed into this function if non-default constants for constructing the uncertainties are desired. This must be done after the parameter sources are properly extracted from the model.\n",
    "\n",
    "## Thermo Uncertainty\n",
    "\n",
    "Each species is assigned a uniform uncertainty distribution in free energy:\n",
    "\n",
    "$G \\in [G_{min},G_{max}]$\n",
    "\n",
    "$dG = (G_{max} - G_{min})/2$\n",
    "\n",
    "Several parameters are used to formulate $dG$.  These are $dG_{library}$, $dG_{QM}$, $dG_{GAV}$, and $dG_{group}$.\n",
    "        \n",
    "$dG =  \\delta_{library} dG_{library} + \\delta_{QM} dG_{QM} +\n",
    "\\delta_{GAV} dG_{GAV} +\n",
    "\\sum_{group} w_{group} dG_{group}$\n",
    "\n",
    "where $\\delta$ is a dirac delta function which equals one if the species thermochemistry parameter contains the particular source type and $w_{group}$ is the weight of the thermo group used to construct the  species thermochemistry in the group additivity method.\n",
    "\n",
    "## Kinetics Uncertainty\n",
    "\n",
    "Each reaction is assigned a uniform uncertainty distribution in the overall ln(k), or ln(A):\n",
    "\n",
    "$d \\ln (k) \\in [\\ln(k_{min}),\\ln(k_{max})]$\n",
    "\n",
    "$d\\ln(k) = [\\ln(k_{max})-\\ln(k_{min})]/2$\n",
    "\n",
    "The parameters used to formulate $d \\ln(k)$ are $d\\ln(k_{library})$, $d\\ln(k_{training})$, $d\\ln(k_{pdep})$, $d\\ln(k_{family})$, $d\\ln(k_{non-exact})$, and $d\\ln(k_{rule})$.\n",
    "\n",
    "For library, training, and pdep reactions, the kinetic uncertainty is assigned according to their uncertainty type.  For kinetics estimated using RMG's rate rules, the following formula is used to calculate the uncertainty:\n",
    "\n",
    "$d \\ln (k) = d\\ln(k_{family}) + \\log_{10}(N+1)*dln(k_{non-exact})+\\sum_{rule} w_{rule} d \\ln(k_{rule})$\n",
    "\n",
    "where N is the total number of rate rules used and $w_{rule}$ is the weight of the rate rule used to estimate the kinetics. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uncertainty.loadDatabase()\n",
    "uncertainty.extractSourcesFromModel()\n",
    "uncertainty.assignParameterUncertainties()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first order local uncertainty, or variance $(d\\ln c_i)^2$, for the concentration of species $i$ is defined as:\n",
    "\n",
    "$(d\\ln c_i)^2 = \\sum_j \\left(\\frac{d\\ln c_i}{d\\ln k_j}d\\ln k_j\\right)^2 + \\sum_k \\left(\\frac{d\\ln c_i}{dG_k}dG_k\\right)^2$\n",
    "\n",
    "We have previously performed the sensitivity analysis.  Now we perform the local uncertainty analysis and apply the formula above using the parameter uncertainties and plot the results.  This first analysis considers the parameters to be independent.  In other words, even when multiple species thermochemistries depend on a single thermo group or multiple reaction rate coefficients depend on a particular rate rule, each value is considered independent of each other.  This typically results in a much larger uncertainty value than in reality due to cancellation error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uncertainty.localAnalysis(sensitiveSpecies, correlated=False, number=5, fileformat='.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Show the uncertainty plots\n",
    "for species in sensitiveSpecies:\n",
    "    print '{}: Thermo Uncertainty Contributions'.format(species)\n",
    "    display(Image(filename=os.path.join(uncertainty.outputDirectory,'thermoLocalUncertainty_{}.png'.format(species.toChemkin()))))\n",
    "    \n",
    "    print '{}: Reaction Uncertainty Contributions'.format(species)\n",
    "    display(Image(filename=os.path.join(uncertainty.outputDirectory,'kineticsLocalUncertainty_{}.png'.format(species.toChemkin()))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlated Uncertainty\n",
    "\n",
    "A more accurate picture of the uncertainty in mechanism estimated using groups and rate rules requires accounting of the correlated errors resulting from using the same groups in multiple parameters.  This requires us to track the original sources: the groups and the rate rules, which constitute each parameter.  These errors may cancel in the final uncertainty calculation.  Note, however, that the error stemming from the estimation method itself do not cancel.  \n",
    "\n",
    "For thermochemistry, the error terms described previously are $dG_{library}$, $dG_{QM}$, $dG_{GAV}$, and $dG_{group}$.  Of these, $dG_{GAV}$ is an uncorrelated independent residual error, whereas the other terms are correlated.   Noting this distinction, we can re-categorize and index these two types of parameters in terms of correlated sources $dG_{corr,y}$ and uncorrelated sources $dG_{res,z}$.\n",
    "\n",
    "For kinetics, the error terms described perviously are $d\\ln(k_{library})$, $d\\ln(k_{training})$, $d\\ln(k_{pdep})$, $d\\ln(k_{family})$, $d\\ln(k_{non-exact})$, and $d\\ln(k_{rule})$.  Of these, $d\\ln(k_{family})$, $d\\ln(k_{non-exact})$ are uncorrelated independent error terms resulting from the method of estimation.  Again, we re-categorize the correlated versus non-correlated sources as $d\\ln k_{corr,v}$ and $d\\ln k_{res,w}$, respectively.  \n",
    "\n",
    "The first order local uncertainty, or variance $(d\\ln c_{corr,i})^2$, for the concentration of species $i$ becomes:\n",
    "\n",
    "\n",
    "$(d\\ln c_{corr,i})^2 = \\sum_v \\left(\\frac{d\\ln c_i}{d\\ln k_{corr,v}}d\\ln k_{corr,v}\\right)^2 + \\sum_w \\left(\\frac{d\\ln c_i}{d\\ln k_{res,w}}d\\ln k_{res,w}\\right)^2 + \\sum_y \\left(\\frac{d\\ln c_i}{dG_{corr,y}}dG_{corr,y}\\right)^2 + \\sum_z \\left(\\frac{d\\ln c_i}{dG_{res,z}}dG_{res,z}\\right)^2$\n",
    "\n",
    "where the differential terms can be computed as:\n",
    "\n",
    "$\\frac{d\\ln c_i}{d\\ln k_{corr,v}} = \\sum_j \\frac{d\\ln c_i}{d\\ln k_j}\\frac{d\\ln k_j}{d\\ln k_{corr,v}}$\n",
    "\n",
    "$\\frac{d\\ln c_i}{d G_{corr,y}} = \\sum_k \\frac{d\\ln c_i}{dG_k}\\frac{dG_k}{dG_{corr,y}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "uncertainty.assignParameterUncertainties(correlated=True)\n",
    "uncertainty.localAnalysis(sensitiveSpecies, correlated=True, number=10, fileformat='.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Show the uncertainty plots\n",
    "for species in sensitiveSpecies:\n",
    "    print '{}: Thermo Uncertainty Contributions'.format(species)\n",
    "    display(Image(filename=os.path.join(uncertainty.outputDirectory,'thermoLocalUncertainty_{}.png'.format(species.toChemkin()))))\n",
    "    \n",
    "    print '{}: Reaction Uncertainty Contributions'.format(species)\n",
    "    display(Image(filename=os.path.join(uncertainty.outputDirectory,'kineticsLocalUncertainty_{}.png'.format(species.toChemkin()))))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
